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Abstract 

In the first part of this paper we present a review of our results concerning the weakly nonlinear 
regime of the mirror instability in the framework of an asymptotic model. This model belongs to 
the class of gradient type systems for which the free energy can only decrease in time. It reveals 
a behavior typical for subcritical bifurcations: below the mirror instability threshold, all localized 
stationary structures are unstable, while above threshold, the system displays a blow-up behavior. It 
is shown that taking the electrons into account (non-zero temperature) does not change the structure 
of the asymptotic model. For bi-Maxwellian distributions functions for both electrons and ions, the 
model predicts the formation of magnetic holes. The second part of this paper contains original 
results concerning two-dimensional steady mirror structures which can form in the saturated regime. 
In particular, based on Grad-Shafranov-like equations, a gyrotropic plasma, where the pressures 
in the static regime are only functions of the amplitude of the local magnetic field, is shown to 
be amenable to a variational principle with a free energy density given by the parallel tension. 
This approach is used to demonstrate that small-amplitude static holes constructed slightly below 
the mirror instability threshold identify with lump solitons of KPII equation and turn out to be 
unstable. It is also shown that regularizing effects such as finite Larmor radius corrections cannot 
be ignored in the description of large-amplitude mirror structures. Using the gradient method, 
which is based on a variational principle for anisotropic MHD taking into account ion finite Larmor 
radius effects, we found both one-dimensional magnetic structures in the form of stripes and two- 
dimensional bubbles when the magnetic field component transverse to the plane is increased. These 
structures realize minimum of the free energy. 
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I. INTRODUCTION 


Magnetic structures in the form of holes or humps associated with maxima or minima of 


plasma density and pressure are often encountered in planetary magnetos’ 


leafs close to both 
3| ) as well. These 


the bow-shock and the magnetopause, and in the solar wind (see e.g. [1 
structures are often viewed as ultra-low frequency (ULF) waves resulting from the mirror 
instability (MI) [41, and, by this reason, called mirror structures. This instability develops 
in a collisionless plasma characterized by a relatively large /3 (a f ew units) and a transverse 
(usually ionic) temperature T± larger than the parallel one Tji, such that the condition for 
mirror instability 

TJT\\ - 1 > /3I 1 (1) 

is fulfilled. Here (3± = 8ttp±/B 2 (similarly, [5 y = 8ttp\\/B 2 ), where p± = nTj_and p\\ = nT\\ 
are the perpendicular and parallel plasma pressures respectively. 

In the Earth magnetosheath, a typical depth of magnetic holes is about 20% of the mean 
magnetic held value and can sometimes achieve 50 %. The characteristic width of such 
structures is of the order of a few ion Larmor radii, and they display an aspect ratio of 
about 7-10. In solar wind, according to 3'], the size of holes may be very different, varying 
from 10 up to 1000 ion gyroradii. In magnetosheath, holes and humps have comparable size, 
and amplitudes. Humps are often observed near the magnetopause where conditions (jT]) 


for development of the MI can be met under the effect of the 
structures are also observed when the plasma is linearly stable [( 


asma compression. Mirror 
7], which may be viewed as 


the signature of a bistability regime resulting from a subcritical bifurcation, whose existence 
was interpreted on the basis of a simple energetic argument within the simplified description 
of anisotropic magnetohydrodynamics j§|. 

The linear mirror instability has been extensively studied both analytically (see es 


lOj), and by means of particle-in-cell (PIC) simulations ll|. As shown in |9j, 12, |13|, the 




instability is arrested at large k due to finite ion Larmor radius (FLR) effects. It turns out 
that wave-particle resonance plays a central role in driving the instability, while the FLR 
effects are at the origin of the quenching of the instability at small scales. In contrast, 
a few years ago, a theoretical understanding of the nonlinear phase remained limited to 
phenomenological modeling of particle trapping jl4, 15] that hardly reproduce simulations 
of Vlasov-Maxwell equations [16 ]. 
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The first nonlinear theory was formulated in 171, [18] where we developed a weakly nonlin¬ 
ear approach to the mirror instability based on the mixed hydrodynamic-kinetic description. 
For the sake of simplicity, an electron-proton plasma with cold electrons was considered 
first. It includes the force-balance equation within the anisotropic MHD and the drift ki¬ 
netic equation for the ions. Close to threshold, the unstable modes have wavevectors almost 
perpendicular to the ambient magnetic field B (k z /k± -C 1) with k±pi -C 1 (pi is the ion 
Larmor radius), so that the perturbations can be described using a long-wave approxima¬ 
tion. The latter allows one to apply the drift kinetic equation (see, e.g. 19|, 20]) to estimate 
the main nonlinear effects that correspond to a local shift of the instability threshold dU). 
All other nonlinearities connected, for example, with ion inertia are smaller. As the result, 
it is possible to derive an asymptotic equation with quadratic nonlinearity of generalized 
gradient type jl7, 18|. The latter property implies an irreversible character of the mirror 
modes behavior, associated with ion Landau damping, where the free energy (or Lyapunov 
functional) can only decrease in time. In this framework, above threshold, the mirror modes 
have a blow-up behavior with a possible saturation at an amplitude level comparable to that 
of the ambient field. Below threshold, all stationary (localized) structures were predicted 
to be unstable. Thus, the system near the MI threshold displays a behavior typical of a 
subcritical bifurcation when the small-amplitude stationary solutions below threshold turn 
out to be unstable; above threshold, the amplitude of magnetic field perturbations tends to 
blow up. It is worth noting that this approach contrasts with the quasi-linear theory 22 1 


that also assumes vicinity of the instability threshold but, being based on a random phase 
approximation, cannot predict the appearance of coherent structures. Phenomenological 
models based on the cooling of trapped particles were proposed to interpret the existence of 
deep magnetic holes jl4, 231These models do not however address the initial value problem 
in the mirror unstable regime. 

i n 

The asymptotic model [171 . [18] was first derived under the assumption of cold electrons. 


Therefore, in our further papers 


24 


25j , we considered how hot electrons can be incorporated 


into the model. The approach we developed is based on the assumption of an adiabatically 
slow dynamics of the mirror structures that allows one to compute the coupling coefficient 
in the weakly nonlinear regime as well as to simplify all calculations of the linear growth rate 
in the case of bi-Maxwellian distributions for both the ions and the electrons. The adiabatic 
hypothesis can be proved perturbatively, and is in particular valid within the asymptotic 
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model. Because this model predicts the existence of subcritcal bifurcation with a blow-up 
behavior above threshold, consistent with the formation of mirror structures with amplitude 
of the magnetic held perturbation comparable with the ambient held, our next step was to 
investigate the properties of possible stationary mirror structures. 

The aim of the present paper is twofolds. The hrst part provides a review of our pre¬ 
vious results concerning the weakly nonlinear model for both cold (13, [l8|) and hot ( 24|) 
electrons. Another goal of this paper is to study steady mirror structures resulting from the 
balance of magnetic and (both parallel and perpendicular) thermal pressures, whose sim¬ 


plest description is provided by anisotropic MHD. 
governed by the Grad-Shafranov (GS) equation 


26 


sotropic MHD equilibria are classically 


271 . |29| . We here revisit this approach 


in the case of anisotropic electron and ion huids where the perpendicular and parallel pres¬ 
sures are given by equations of state appropriate for the static character of the solutions. 
However, the MHD stationary equations, at least in the two-dimensional geometry, turn out 
to be ill-posed. As a consequence, these equations require some regularization. As done 
in a similar context of pattern formation 30], an additional linear term involving a square 
Laplacian is added. For nonlinear mirror modes, regularization can originate from finite 
Larmor radius (FLR) corrections, which are not retained in the present analysis based on 
the drift kinetic equation (see, e.g. {l7, 1^|). 

The paper is organized as follows. In Section II, we discuss the linear mirror instability 
near the MI threshold. Section III is devoted to the derivation of weakly nonlinear asymp¬ 
totic model, in the simplest case of cold electrons, and to its properties, including possible 
stationary states (below the MI threshold) and blow-up behavior (above threshold). Section 
IV deals with accounting electrons in the asymptotic model. Here we develop the adiabatic 
approach for finding contributions from electrons to both the linear growth rate and the 
nonlinear coupling coefficient entering the asymptotic model. In Section V, we formulate 
the variational principle for the stationary anisotropic MHD when both parallel and trans¬ 
verse pressures depend on the magnetic field amplitude with a free energy given by the space 


integral of the parallel tension. In this case, as well known 


31 


-36], the parallel component 


of the MHD equation is satisfied identically. In Section 6, the anisotropic Grad-Shafranov 
equations are revisited when the gyrotropic pressures depend only on the local magnetic held 
amplitude that, as shown in the forthcoming sections, is specific of nonlinear mirror modes. 
In this case the stationary anisotropic MHD represents an hydrodynamic integrable-type 
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system and for this reason requires the renormalization due to FLR effects. In this Section, 
it is shown also that the equations of state resulting from an adiabatic approximation of 
the drift kinetic description, require a regularization because of an overestimate of the con¬ 
tributions from the particles with a large magnetic moment. We discuss in particular the 
small-amplitude regime and show that the pressure-balanced structures are then governed 
by the KPII equation which possesses lump solutions. Numerical simulations reproduce 
these special structures, that turn out to be unstable. Computation of stable solutions lead 
to large-amplitude purely one-dimensional solutions in the form of stripes that appear to be 
sensitive to the regularization process, an indication that the regime cannot be captured by 
the drift kinetic approximation and that finite Larmor corrections and trapped particles are 
to be retained. Section VII aims for presentation of the numerical results for two-dimensional 
(depending on x and y coordinates) stationary mirror structures when the magnetic held 
B has also a B z component. In particular, we show that for small B z stationary structures 
realizing the minimum of the free energy, below and above the threshold, have the form of 
stripes which are one-dimensional structures with constant magnetic held outside and inside 
the stripes. The transient region, between outer and inner regions, for the stripes represents 
the magnetic well which structure is defined by the FLR contributions to the free energy. 
With increasing B z , instead of stripes, the free energy has its minimum for bubble-type 
structures with an elliptic form. When B xy —>■ 0 these bubbles become circular. In this 
case, FLR effects play a role of the surface tension. Section VIII is the conclusion. 

II. MAIN EQUATIONS AND MIRROR INSTABILITY 

Consider for the sake of simplicity, a plasma with cold electrons. To describe the mirror 
instability in the long-wave limit it is enough to use the drift kinetic equation for ions 
ignoring parallel electric held E y and transverse electric drift: 

^ +B|b .V/-Mb-VB^- = 0. (2) 

In this approximation ions move along the magnetic held (b = B /B) due to the magnetic 
force p b ■ Vi? where p = vj_/(2B) is the adiabatic invariant which plays the role of a 
parameter in equation (J2]) . Both pressures p« and p± are given by 

p\\ = mB / v%fdpdv\\dip = m vi j/d 3 u, (3) 
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Equation (J2]) with relations ([3]) and (j4|) are supplemented with the equation expressing the 
balance of forces in a plane transverse to the local magnetic held 


n{-v( P x + fd 


47T 


+1 1 + bskp-l - PI 


B VB 

47T } 


1 = 0. 


(5) 


Here, consistently with the long-wave approximation, we neglect both the plasma inertia 
and the non-gyrotropic contributions to the pressure tensor. Furthermore, H^ = So- — bibk 
denotes the projection operator in the plane transverse to the local magnetic held. In this 
equation, the hrst term describes the action of the magnetic and perpendicular pressures, 
the second term being responsible for magnetic lines elasticity. 

The equation governing the mirror dynamics is then obtained perturbatively by expanding 
Eqs. ([2]) and d5]) . In this approach, the ion pressure tensor elements are computed from the 
system (J2]), (J5]), near a bi-Maxwellian equilibrium state characterized by temperatures T± 
and Tji and a constant ambient magnetic held Bo taken along the ^-direction. 

From Eq. ([5]) linearized about the background held B 0 by writing B = B n + B (B 0 B) 


with B ~ e ia,t+lk ' r ; we have 


P { ± 


BqBz 


47r 


K 

Afj 


1 + 


/3_l — f 3 1| \ B 0 B Z 


47r 


( 6 ) 


Here k z and k± are the projections of the wave vector k, and is calculated from the 
linearized drift kinetic equation (J2|): 

d fW dfM dB z dfW _ . 
dt + ' 11 dz 11 dz dv\\ 

In Fourier space, this equation has the solution 


/ (1, = 


pB z 


k 


6 >/(°) 


uj — k z v || " dv 


(7) 


The mirror instability is such that uj/k z <C v t h\\ = ^Ty/m. This means that the ions con¬ 
tributing to the resonance c o — k z V\\ = 0, correspond to the maximum of the ion distribution 
function. 














After substituting (J7J) into the first order term for perpendicular pressure 
forming integration, we get 




0±\ B 0 B Z iy/nu /3 2 ± B 0 B Z 


and per- 


( 8 ) 


4vr \k z \v t h\\ /3|| An 

The hrst term in (J 8 ]) is due to the difference between perpendicular and parallel pressures, 
while the second one accounts for the Landau pole. 

Equation (jHJ) together with (J2J) yield the growth rate for the mirror instability in the drift 
approximation where FLR corrections are neglected 4] 


7 = \kz\v th \\ 




V^/Sj 


P ±_ 1 _ J_ 


hi 


-x 


(9) 


,0|| ‘ 0u k 2 ± /3 ± 

where y = 1 + (0j_ — 0 1 |)/ 2 . The instability takes place when the criterion (0Q) is fulfilled and, 
near threshold, develops in quasi-perpendicular directions, making the parallel magnetic 
perturbation dominant. 

As shown in Refs. 9, 12, |l3], when the FLR corrections are relevant, the growth rate is 
modified into 


7 = \k z \vth\\ 

where £ = 0 _lX _1 (0_l/0|| - 1 


0IA 


M. 

k 2 


-kip 2 


( 10 ) 


\Ar0i L k'i Ax' 

0J 1 ) and the ion Larmor radius p z = v t h±/co C i is defined with 
the transverse thermal velocity Vth± = \J^T^/m and the ion gyrofrequency co ci = eB 0 /(mc). 
This growth rate can be recovered by expanding the general expression given in |9|. in the 
limit of small transverse wavenumbers, ft can also be obtained directly from the Vlasov- 
Maxwell (VM) equations in a long-wave limit which retains non gyrotropic contributions 0 - 
ft is important to note that the expression (TTUT) for 7 is consistent with the applicability 
condition ou/k z v t h\\, he. when the supercritical parameter |e| <C 1. In this case the 
instability saturation happens at small k± oc y/e due to FLR and for almost perpendicular 
direction in a small cone of angles, k z / k± oc y/e. As a result, the growth rate 7 oc £ 2 , so 
that, when defining new stretched variables by 


k z = eK z pi 1 (2/73)y 1/2 , 

k± = y/s(2/V3)K ± pr 1 x 1/2 , ( 11 ) 

7 = £ 2 T(2/73)12 ( 0 F 0 ± )~ 1 (x 0 |,/ 0 ±) 3/2 , 

it takes the form 

T = \K Z \ (l - K 2 z fKl - Kl) . (12) 
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Hence it is seen that, in the (K± — 0) plane (0 = K z /K±), the instability takes place inside 
the unit circle: @ 2 + K 2 < 1. The maximum of T is obtained for K± = 1/2, @ = ±1/2 and 
is equal to T max = 1/8. Outside the circle, the growth rate becomes negative (in agreement 


with 


id)- 


III. WEAKLY NONLINEAR REGIME: ASYMPTOTIC MODEL FOR COLD 

ELECTRONS 


A. Derivation 


As it follows from (15|) , in the linear regime, near the instability threshold, the fluctuations 
of perpendicular and magnetic pressures almost compensate each other (compare with (J9j)). 
Therefore, in the nonlinear stage of this instability, we can expect that the main nonlin¬ 
ear contributions come from the second order corrections to the total (perpendicular plus 
magnetic) pressure, i.e. 


P± + 


BnB- 


Bz 


in + Pl + 8? “ _X A 1 in 
This result can be obtained rigorously by means of a multi-scale expansion based on the 
linear theory scalings (1111) . For this purpose, we introduce a slow time T and slow coordinates 
R in a way consistent with (HTTh and expand the magnetic field fluctuations as a powers 
series in e 1 / 2 : 


<9? B n B~ 


(13) 


B z = eB™ + 0{e 2 ), Bi = £ 3 / 2 b|l 3/2) + 0(e 5/2 ), 


(14) 


where B^ 2 ) are assumed to be functions of R and T. Using these expressions, it is easy to 
establish that quadratic nonlinear terms coming from the expansion of n in (J5|) as well as 
from the second term in the r.h.s. of Eq. (1T31) are small in comparison with the quadratic 
term originating from the magnetic pressure in Eq. (1T31) . Thus, to get a nonlinear model 
for mirror dynamics, it is enough to End p'^. The expansion (TT4l) induces a corresponding 
expansion for the distribution function and for both pressures. Defining 


= nm J v±f^vj_dv±dv 


from ([4j) we have 

P± ] = {B^/B.fpf + 2 (BW/Bo)^ +f%\ 
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up to an additional contribution proportional to that cancels out in the final equation 
due to the threshold condition. 

On the considered time scale, the effect of nonlinear Landau resonance is negligible in 
the contribution to f ( - 2> that can thus be estimated from the equation 


V\\ 


<9/( 2 ) 

dz 


(VA&i )ei 1) 


5/(0) 
dz dv\\ 


0 . 


For an equilibrium bi-Maxwellian distribution, we have 


/ (2) = (2 M 2 /4 || )(B«) 2 /<»> (15) 

and thus 

pf = (p ± - 4131/P, + 3 01/$1) % 

As a consequence, because of the vicinity to threshold we obtain 

# + g-(l + #)g>0. ( 16 ) 

Then rewriting equation (1T3T) using the slow variables (HTT) and rescaling the amplitude 


B z /B 0 = e2xP±(l + Ai) u i 


we arrive at the equation H Q 

du 


dT 


= Kp 


a ~ A± 'jzi + A± ) “ _ 3 “ 2 


( 17 ) 


Here a = ±1, depending of the positive or negative sign of e, Iiz = —Bdz is a positive 
definite operator (whose Fourier transform is \Kz\) 1 H is Hilbert transform: 


Hf{Z) = -VP 


7T 


f°° f(Z') 

Too Z'-Z' 


dZ'. 


As seen from the equation, its linear part reproduces the growth rate (fT2]h In particular, 
the third term in the r.li.s. accounts for the FLR effect. 

Equation (fT71) simplihes when the spatial variations are limited to a direction making a 
fixed angle with the ambient magnetic field. After a simple rescaling, one gets 


du 

dT 


= K- 


a + mp) u - 3u2 


(18) 


where S is the coordinate along the direction of variation. This equation can be referred to 
as a “dissipative Korteveg-de Vries (KdV) equation”, since its stationary solutions coincide 
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with those of the usual KdV equation. The presence of the Hilbert transform in Eq. (1T81) 
nevertheless leads to a dynamics significantly different from that described by soliton equa¬ 
tions. Besides, it is worth noting also that Eq. (1171) in the two-dimensional case has some 
similarity with the KP equation (see, Section IV). 


B. Properties of the asymptotic model 


Equation (H71) (and its ID reduction (USD as well) possesses the remarkable property of 
being of the form 

du _ 6F 


where 


a 2 , 1„. A-lo2„. , 1 /V7 „.\2 , „.3 


'-u 2 + ^uA^d\u + ^ (V ± m) z + u A 


F = J 

= —<j N/ 2 + Ji/2 + /2/2 + / 3 


dR 


(19) 


has the meaning of a free energy or a Lyapunov functional. This quantity can only decrease 
in time, since 


-dR = — 


SF-^ 5F 


Kz—dFi, < 0 . 


( 20 ) 


dF _ fSFBu 

dt J 5u dt J 5u 5u 

This derivative can only vanish at the stationary localized solutions, defined by the equation 


^ = (a - AT 1 ^ + A ± ) u - 3u 2 = 0. 


5u 


dZ 2 


( 21 ) 


We now show that non-zero solutions of this equation do not exist above threshold (a = 
+1). For this aim, following Ref. js8|, we establish relations between the integrals N , R, I 2 
and I 3 , using the fact that solutions of Eq. (12T1) are stationary points of the functional F 
(i.e. SF = 0). Multiplying Eq. (12T1) by U and integrating over R gives the first relation 


(j N — I\ — I 2 — 3/3 — 0. 


Two other relations can be found if one makes the scaling transformations, Z —>■ aZ, R^ —>■ 
bH±, under which the free energy (TT9l) becomes a function of two scaling parameters a and 
b 

F(a, b) = —+ ~^ a + 7 3 afe 2 . 
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Due to the condition 8F = 0, the first derivatives of F at a = b = 1 have to vanish: 


dF aN h J 2 r 

&=--- 2 + 2 +/3 
dF 

— = - a N + 2/i + 2/3 


0, 

0. 


Hence, after simple algebra, one gets the three relations 

h + — N = 0, I 3 = —2/i, J 2 = 3/i. 


For a — +1, the first relation can be satisfied only by the trivial solution u = 0, because 
both integrals I\ and N are positive definite. In other words, above threshold, nontrivial 
stationary solutions obeying the prescribed scalings do not exist. 

In contrast, below threshold, stationary localized solutions can exist. For these solutions, 
the free energy is positive and reduces to F s = N/2. Furthermore, J 3 = / U 3 d 3 R < 0 which 
means that the structures have the form of magnetic holes. As stationary points of the 
functional F, these solutions represent saddle points, since the corresponding determinant 
of second derivatives of F with respect to scaling parameters taken at these solutions is 
negative (d aa FdbbF — ( d a bF ) 2 = —2N 2 < 0). As a consequence, there exist directions in the 
eigenfunction space for which the free-energy perturbation is strictly negative, corresponding 
to linear instability of the associated stationary structure. This is one of the properties for 
subcritical bifurcations. 

As a consequence, starting from general initial conditions, the derivative dF/dt (1201) is 
almost always negative , except for unstable stationary points (zero measure) below thresh¬ 
old. In the nonlinear regime, negativeness of this derivative implies / u 3 d 3 R < 0, which 
corresponds to the formation of magnetic holes. Moreover, this nonlinear term (in F) is 
responsible for collapse, i.e. formation of singularity in a finite time. 


C. Blow-up 


In order to characterize the nature of the singularity of Eq. (1151) . it is convenient to 
introduce the similarity variables £ = (To — T)~ 1 // 3 S, r = — log(To — T), and to look for a 
solution in the form U = (T 0 — T)~ 2 ^ g(£, r), where g(£,r) satisfies the equation 


99 2 , £ dg _ — 
fr + 3 9+ 3dt- R t 


&g_ 

8 ^ 


3y 


+ e~ r iqg. 
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As time T approaches To (r —> oo), the last term in this equation becomes negligibly small 
and simultaneously d T g —» 0 so that asymptotically the equation transforms into 


3 9 


£dg p 

aH A{ 


d 2 a o 2 


( 22 ) 


For the free energy this means that close to To the first term ~ N turns out to be much 
smaller in comparison with all other contributions, in particular with / f/ 3 dS. 

At large |£|, that corresponds to the limit T —* T 0 , the asymptotic solution g of Eq. (f22ll 
obeys 

29+f = cr 2 

where C = | 5 ,2 (£ , ) C ^Z > 0; and has the form g— C£,~ 2 log |£/£ 0 |- For U, it gives the 

asymptotic solution 

C 

Uasymp 1—2 lo§ I'Z 1 —‘0 (^) | 

with 5 0 (t) = (T 0 — T) 1 / 3 ^, that, as T —>■ To, has an almost time independent tail. For 
|S| < (To — T ) 1 ^ 3 |£o|, the solution is negative and becomes singular as 5 approaches the 
origin. 

Asymptotically self-similar solutions can also be constructed in three dimensions, when 
rescaling the longitudinal coordinate by (T 0 — T) 1 / 2 , the transverse ones by (T 0 — T ) 1 / 4 and 
the amplitude of the solution by (To — T)” 1 / 2 . Existence of a finite time singularity for the 
initial value problem can be established for initial conditions for which the functional F is 
negative, when the term involving a can be neglected, an approximation consistent with the 
dynamics: 

h h 

F —> Tii m = — + — + / 3 . (23) 

To prove this statement, consider the operator K~ l , (inverse of the operator K z ), which is 
dehned on functions obeying / U(Z, Rj _)dZ = 0, a condition consistent with Eq. (ITTlh Then 
the time derivative of Fy un can be rewritten through the operator JCT 1 as follows, 


= ~ f U T K^U T dR < 0. (24) 

Consider now the positive definite quantity N = J UK^UdR > 0, whose dynamics is 
determined by the equation 


dN 

Hr 


—2 (Ii + I 2 + 3 / 3 ) — — 6 Fy m + Iy + I 2 . 


(25) 
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Let F] im be negative initially, then at T > 0 the r.li.s. of (1251) will be positive, and, as a 
consequence, N will be a growing function of time. 

Introduce now the new quantity S = —Fu m /N which is positive definite if F Hm \ T=0 < 0. 
The time derivative of S is then defined by means of Eqs. (1231) and (1231) : 


dS 

dT 


F]\ m N r 1 
N 2 + IV 


UrK^UrdR. 


(26) 


The second term in the r.h.s. of this equation can be estimated using the Cauchy- 
Bunyakowsky inequality: 

UK^UrdR < 2 N 1/2 ^ j U T K; x U T dR 




that gives 


U T K- l U T dR > N^/{AN). 


Substituting the obtained estimate into Eq. (1261) and taking into account definition (1231) for 
F\ im and Eq. (1251) as well, we arrive at the differential inequality for S (compare with |39|): 


dS_ > Nt 
dT ~ N 2 


n t 

— ~ Fim 


> 15 S 2 . 


Integrating this first-order differential inequality yields 


S > 


1 

15 (T„-T)' 


(27) 


Here the collapse time To = (ISS'o)” 1 is expressed in terms of the initial value S\t=o = So. 
It is interesting to mention that the time behavior of S given by the estimate (1271) coincides 
with that given by the self-similar asymptotics. 


D. Conclusion of Section III 

We have presented an asymptotic description of the nonlinear dynamics of mirror modes 
near the instability threshold. Below threshold, we have demonstrated the existence of 
unstable stationary solutions. Differently, above threshold, no stationary solution consistent 
with the prescribed small-amplitude, long-wavelength scaling can exist. For small-amplitude 
initial conditions, the time evolution predicted by the asymptotic equation (TT7T) leads to a 
finite-time singularity. These properties are based on the fact that this equation belongs 
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to the generalized gradient systems for which it is possible to introduce a free energy or a 
Lyapunov functional that decreases in time. 

The singularity formation as well as the existence of unstable stationary structures below 
the mirror instability threshold obtained with the asymptotic model, can be viewed as 
features of a subcritical bifurcation towards a large-amplitude state that cannot be described 
in the framework of the present analysis. Such an evolution should indeed involve saturation 
mechanisms that become relevant when the perturbation amplitudes become comparable 
with the ambient held. 


IV. ADIABATIC APPROACH: ACCOUNT OF ELECTRONS 


The mirror instability, as known, is a kinetic instability whose growth rate was first ob¬ 
tained under the assumption of cold electrons j4|, a regime where the contributions of the 
parallel electric held En can be neglected. However, in realistic space plasmas, the electron 


temperature can hardly be ignored 


42]. The linear theory retaining the electron temperature 


and its possible anisotropy, in the quasi-hydrodynamic limit (which neglects finite Larmor 
radius corrections), was developed in the case of bi-Maxwellian distribution functions by 
several authors (see e.g. [43 l |, 9], 10|). A general estimate of the growth rate under the sole 
condition that it is small compared with the ion gyrofrequency (a condition reflecting close 
vicinity to threshold) is presented in 2^]. Like for the cold electrons case, the instability 
develops in quasi-perpendicular directions, making the parallel magnetic perturbation dom¬ 
inant. This analysis includes in particular regimes with a significant electron temperature 
anisotropy for which the instability extends beyond the ion Larmor radius. In the limit 
where the instability is limited to scales large compared with the ion Larmor radius pi, only 
the leading order contribution in terms of the small parameter ^/ (\k\ z v\\i) is to be retained 
in estimating Landau damping, and the growth rate is given by 
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Here, Tj_ a and Tju are the perpendicular and parallel (relative to the ambient magnetic held 
B 0 taken in the z direction) temperatures of the species a (a — i for ions and a = e for 
electrons ), 9± = Tj_ e /T ±i , 9\\ = Tj| e /Tj|j and /3± = /3±i+/3± e with /3j_ a = 8np± a /B^ where p± a 
is the perpendicular thermal pressure (similar definition for /3||). Furthermore, the parallel 
thermal velocity is defined as v\\ a = <j2T\\ a /m a , and p* = (2Tj_i/m p ) 1 / 2 /Qi denotes the ion 
Larmor radius (hh = eBo/rriiC is the ion gyrofrequency). 

The growth rate given by Eq. (1281) has the same structure as in the cold electron regime 


considered in the previous sections, and given first time in 13l in the case of bi-Maxwellian 

n n 

ions and then generalized in [9] and [10] to an arbitrary distribution function. The first term 
within the curly brackets provides the threshold condition which coincides with that given 
in J42|,jl3 1, 441. The second one reflects the magnetic held line elasticity and the third one 


(where F depends on the electron temperatures due to the coupling between the species 
induced by the parallel electric held which is relevant for hot electrons) provides the arrest 
of the instability at small scales by finite Larmor radius (FLR) effects. 


A. Asymptotic model for hot electrons 

Now we extend to hot electrons the weakly nonlinear analysis developed for cold electrons 
in the previous section. Since in this asymptotics, FLR contributions appear only at the 
linear level, the idea is to use the drift kinetic formalism to calculate the nonlinear terms. 
We show that the equation governing the evolution of weakly nonlinear mirror modes has 
the same form as in the case of cold electrons. In particular, the sign of the nonlinear 
coupling coefficient that prescribes the shape of mirror structures, is not changed, in the case 
of bi-Maxwellian distributions for both electrons and ions, but can be changed for another 
distributions. This equation is of gradient type with a free energy (or a Lyapunov functiona’ 


which is unbounded from below. This leads to finite-time blowing-up solutions 
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associated with the existence of a subcritical bifurcation 17], |l8j. To describe subcritical 
stationary mirror structures in the strongly nonlinear regime, we present an anisotropic 
MHD model where the perpendicular and parallel pressures are determined from the drift 
kinetic equations in the adiabatic approximation, in the form of prescribed functions of the 
magnetic field amplitude only. 

A main condition governing the nonlinear behavior of mirror modes is provided by the 
force balance equation 

(B • V)B 


( IS- 1 


47T 

+ 80 
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1 + ^{p±-p\\) 


47T 


+B 


(B • V) ( 


P± 


V B 2 


v-n = o, 


(30) 


where a gyroviscous contribution II originating from FLR effects (compare with (11711). Note 
that FLR contributions also enter the gyrotropic pressures. Here the pressure tensor and its 
components are viewed as the sum of the contributions of the various species. In particular 
p± = J2aP±a and p\\ = J2aP\\a- When concentrating on scales large compared with the 
electron Larmor radius, the non-gyrotropic correction II to the pressure tensor originates 
only from the ions. As mentioned above, it is enough to retain this contribution only at 
the linear level with respect to the amplitude of the perturbations. As in the case of cold 
electrons, the other linear and nonlinear contributions can be evaluated from the drift kinetic 
equation 

^ + ^b.V/„ + (-Mb.VB + ^)^ = 0 (31) 

for each type of particles. 

We ignore the transverse electric drift which is subdominant for mirror modes. In this 
approximation, both ions and electrons move in the direction of the magnetic field under the 
effect of the magnetic force p b • VB and the parallel electric field E\\ = —b ■ V0 where the 
magnetic moment p = v\/(2B) is an adiabatic invariant which plays the role of a parameter 
in Eq. (13T|) . Here 0 is the electric potential. The quasi-neutrality condition n e — n t = n, 
where n a — B f f a dpdv\\dp = / f a d 3 v, is used to close the system and eliminate E\\. 

In this framework where FLR corrections are neglected, the gyrotropic pressures p\\ a and 
pj_ a are given in terms of the corresponding distribution functions f a by 


Pa || = m a B j 
Pa .l = m a B 2 


f\fadpdv\\d(p 



pf a dpdv\\d(p = -7Ti, 


V\\fad 3 V, 

J v 2 ± f a d 3 v. 
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The equation governing the mirror dynamics is obtained perturbatively by expanding 
Eqs. (I30|) . (l3TTi and the quasi-neutrality condition. In this approach, the pressure tensor 
elements for each species are computed near a bi-Maxwellian equilibrium state characterized 
by temperatures Tj_ a and T\\ a and a constant ambient magnetic field B 0 taken along the 
^-direction. 

B. Linear instability 

Before turning to the nonlinear regime, we briefly reformulate the linear theory in the 
framework of the drift kinetic approximation, in order to specify the notations. 

From Eq. (f30l) , linearized about the background field B 0 by writing B = B 0 + B (B 0 
B) with B ~ e-iut+i kr^ we arr j ve (|5|) . where p"^ has to be calculated from the 

linearized drift kinetic equation (T3TT) after elimination of the parallel electric field using the 
quasi-neutrality condition. Note that as for the case of cold electrons, near the instability 
threshold the leading terms in (J6]) corresponding to perturbations of perpendicular and 
magnetic pressures are compensated by each other and therefore one needs to retain the 
next order terms responsible for both elasticity of magnetic field lines and FLR corrections. 

The linearized drift kinetic equation reads 



(32) 


where we assume each f^ to be a bi-Maxwellian distribution function 



(33) 


with A a = n 0 ‘m a /(27ry/TTV\\ a T ±a ). 

In Fourier representation, Eq. (1321) is solved as 


U! — K Z V\\ OV || 

The neutrality condition in the linear approximation reads 

J f' 1 - 1 dv\\dpd<p = v\\dpd(p, 

that allows one to express the potential 0 in terms of B z . We have 



(34) 



(35) 



(36) 


19 

















Here we assume that oj/k z v\\ z = ^‘2T\\ l /m ll so that the contribution from the Landau 
pole is small (£ = \Jntjj / (\k z \v\\i) <C 1 ). 

An analogous calculation for the electrons, neglecting the contribution of the correspond¬ 
ing Landan resonance because of the small mass ratio (assuming the ratio of the electron to 
ion temperatures is not too small), gives 

B z 


j f 7 dv z dpd<p = 


n 0 

BnT\ 


0-t lie 




The quasi-neutrality condition then reads 
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and leads to the estimate 
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Thus, for cold electrons (6_ j_ — 9\\ = 0), (j) vanishes and the influence of the parallel electric 
held on the mirror instability becomes negligible. Interestingly, when 0j_ = 9 m, only the 
Landan pole contributes to 


T ±i 9\\ B z 

ecp «- -it —. 

^ 1 + 011 Bo 
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we get 


y rn„ J /./('V//. *1|V « -nc^A (c + i(D), 


(42) 


where the coefficients C and H, dehned above, are both positive. In the cold-electron limit, 
C — > 2 and D 2. 

(&± — 0||) / \ (2 + 0_L + ^||) 

It is worth noting that the terms — -+--+— in C and [9± — 0m) —- 7 H- in 

«l(l + »l) V IM (l + 9„) 2 

H originate from the contributions of the electrostatic potential (j) to p 7 , and vanish for 
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0 _l = 6 1 |. Furthermore, in this limit, the real part of the perpendicular pressure fluctuations 
is the sum of two independent contributions originating from the ions and the electrons. 
Differently, only the ion Landau pole contribution is retained in the imaginary part. We 
finally get 


P± ] = ~^n 0 T ±i 
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Substituting this expression into the linearized force balance equation yields 


T 


n Q T 1 _~Dii = 2n Q T u {\ + e A _) 

J- 11 7. 


X 


-± i 


C + — + 


kl 


2Tjij(l + 9±) p ± kl/3 ± 


-X 


and thus the linear growth rate 
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up to the FLR term which is not captured 
by the drift kinetic approximation. As 9 —y 0 , the growth rate reduces to the usual form 
given 


in 4] 
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In the presence of hot electrons, the mirror instability arises when 

T±i 1 


r = 


2Tj|j(l + 6±) Q\\(Q\\ + 1) L 


+ 0i) 2 + 20,|(0i + i) 


(46) 


and, near threshold, develops in quasi-perpendicular directions, making the parallel magnetic 
perturbation dominant. This instability condition can be also rewritten in the form given 


m 


43|. 


Note that the growth rate derived above is valid provided the condition 7 /k~ u t h||i is 
fulfilled. Furthermore, the instability is arrested by FLR effects at scales that are too small 
to be captured by the drift kinetic asymptotics. 
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C. General pressure estimates 


As demonstrated in Section 2 (see also jl7, 181), the scalings (TTll) resulting from the linear 
theory near threshold, when k z and kj_ vary proportionally to e and \fe respectively, while 
the instability growth rate behaves like ~ e 2 , imply an adiabaticity condition, or, in another 
words, this leads to the stationary kinetic equation 


U||b ' V/ a 


(b-V) 


/nB + — 
m a j 


dfg 

dvn 


0 . 


(47) 


It in fact turns out that Eq. (1471) is exactly solvable, the general solution being an 

arbitrary function of all integrals of motion f a = g a (n,W a ,q ) of the particle energy 

e a 

W a = — + fiB H- —(j), of /i and of variables q responsible for labeling the magnetic held 

2 m Q 

lines. As we see in the previous case on the example of cold electrons, the dependence on 
q does not appear in the weakly nonlinear regime, analyzed perturbatively. In the next 
section, we will return to this question and discuss it in more detail. Below we will ignore 
this dependence, considering only the case when f a has two arguments /i and W a . 

To find the function g Q (/u,W a ) in this case, we use the adiabaticity argument which 
means that, to leading order, g a as a function of its arguments /i and W a retains its form 

during the evolution. Therefore, the function g a (fi, W a ) is found by matching with the initial 

v 2 

distribution function given by Eq. (1331) and corresponding to 0 = 0 and W a = -j + fiB 0 . 
We get 
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Thus, W a ) is a Boltzmann distribution function with respect to W a but, at hxed W a , 
it displays an exponential growth relatively to /i if T± a > T\\ a . This effect can however be 
compensated by the dependence of W a in /i. This means that only a fraction of the phase 
space (/i, W a ) is accessible, a property possibly related with the concepts of trapped and 
untrapped particles. 
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Note that expanding Eq. (1481) relatively to B z and (jZB reproduces the first order contri¬ 
bution to the distribution function (1341) with oj = 0 and also the corresponding exp ression 
for the second order correction (fl5]l found in the previous section (see also (l7|, jl8j) in the 
case of cold electrons. It should be emphasized that Eq. (1481) only assumes adiabaticity and 
remains valid for finite perturbations. 

The function g a can also be rewritten in terms of nil, u_l and </> as 

r m„v 2 
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which can be viewed as the bi-Maxwellian distribution function with the renormalized trans¬ 
verse temperature 
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Note the Boltzmann factor exp — [e a (j)/T\\, 
ion distribution function was obtained in 


] in the expression of g a . For cold electrons, the 


51] by assuming that the distribution remains bi- 


Maxwellian and owing to the invariance of the kinetic energy and of the magnetic moment. 
This estimate obtained by neglecting both time dependency (and consequently the Landau 
resonance) and finite Larmor radius corrections reproduces the closure condition given in 

fl. 

After rewriting Eq. (1481) in the form 
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the quasi-neutrality condition gives 
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Interestingly, the electron density 


n e = n 0 
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5n 


1 + 


T\ P B — Bn 


-l 


Tii, 


B n 


exp 


(52) 


has the usual Boltzmann factor exp e0/Xj| e and also an algebraic prefactor depending on 
the magnetic held B. In the case of isotropic electron temperature (T± e = Tj \ e = T e ), the 
electron density has the usual Boltzmann form n e = no exp [ecp/T e ]. 

The above formula for 0 shows that the potential vanishes in two cases: for cold electrons 
and when electron and ion temperature anisotropies a e and a* are equal, a case first time 
mentioned in the linear theory of the mirror instability [12l . ll3L 143]. 

Equation (j5Tj) allows one to evaluate explicitly the perpendicular pressure for each species 


P±a = rn a B 2 J pg a d[idv\\d(p 
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where exp is given by Eq. (ED- 


Hence, simple algebraic procedure gives the following expression for the parallel pressure 

Q. H: 

l 1i 

Pi = "oPlK + ^) (1 + aear . (1+aiMr . . (53) 

where u = B/B 0 — 1, a a = Tj_ a /T|| a is the parameter characterizing the anisotropy of 
distribution function f a , and c a = Tj|Q,(Tj| e + T^)^ 1 in the case of a proton-electron plasma. 
As it will be shown in the next section, the perpendicular pressure can be easily found by 
means of the general relation 
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Substitution of (1531) into this expression yeilds 
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Hence one can see that both pressures have the singularities at u = —a” 1 corresponding to 
the magnetic held 

n — 1 

(56) 


B, = i :„-—^ < /;, 


In the limiting case of cold electrons, p\\ = noT||(l + u)(l + au)^ 1 displays a pole singularity. 
Here, Ty and the anisotropy parameter a correspond to ions only. Such an equation of state 
was previously derived by a quasi-normal closure of the fluid hierarchy 
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The above singularities are presumably related to an overestimated contribution from 
large /i, corresponding either to small B or to large a transverse kinetic energy. In both cases, 
the applicability of the drift approximation breaks down and we are thus led to introduce 
some cut-off type correction near fji* a . In a simple variant, we take f a = C a exp(— m a W a /T\\ a ) 
at /i > /x*, with some positive constant C a , and f a retains its original form (H51) for ji < /x*. 
For cold electrons, the parallel ion pressure is modified into p\\ = tiqT\\G{B, r) with 


G(B,r) 
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exp[~r(B - B a )] -1 
exp[—r(i? 0 - B s )} - 1 

Here, C is a (small) constant, and r = m/x*/T\\. Noticeably, regularization leads to a non¬ 
singular positive pressure for all B , including when B —)■ 0. The modification for p» in the 
case of hot electrons is not specified here because the expressions are algebraically much 
more cumbersome but do not involve any additional difficulty. 

After these remarks, one can easily derive the asymptotic model with account of hot 
electrons. The basic idea is the same as we used already while derivation the model (fTD) for 
cold electrons. To derive the asymptotic model, we can of course forget about renormaliza¬ 
tion of the function G(B, r ) because we need to consider the expansion of with respect 
to small amplitude u by taking into account in this expansion only the second term ~ u 2 
which defines the nonlinear coupling coefficient for (H71) . For (155|) the quadratic contributions 
originating from + (B — Bq) 2 /(8 tt) are collected in a term A ( B gf 13 ) with 



A = 77, 0 |Tj_j^3a^ — 4aj + 1 

+Ci(a e Cli'j ~(1 T Cii ) 2 + 3 di ^ 


+Tj_ e (3a- — 4a e T 1 T c e (o e — Oj) 

(1 + c e )(a e — cii ) + 2 — 3a e ) | + 


g| 

8 tt 


B 2 

The value A c of A at threshold is obtained by expressing —- by means of Eq. 
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where A c = A c /(n 0 T_u). 

Supplementing the corresponding quadratic terms in Eq. 
expansion, to the dynamical equation 
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that extends the result of |l7l . [37] valid for cold electrons. As demonstrated in 
sign of the nonlinear coupling A c defines the type of subcritical structures, namely holes 
(A c > 0) or humps (A c < 0). It turns out that the sign of the nonlinear coupling can be 
determined analytically in a few special cases. 

(i) Limit 9\\ <C 0±: 

A 3 E" (59) 
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(ii) Equal anisotropies (9± — 9\\) 
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(in) Isotropic electron temperature: The coefficient A c can be rewritten in the form 


A c = n 0 (ai - l){Tj_i((3ai - 1) 


+Cj 


-(1 + a) (aj — 1) + 2 — 3a* J 


1 B ^ 

+T e c e — (1 + c e ) (a* — 1) + 1 } + -2-. 

L 2 J 07r 


Furthermore, at threshold 
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Hence, we simultaneously have two inequalities a* > 1 and Tj_ e c e > Tj_j(cj — 2). Therefore, 
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which is positive because 1 


Ci — Cp 
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+}. 

> 0 and ai — 1 > 0. 


(iv) More general conditions: A numerical approach was used. For this purpose it is of 
interest to display in Fig. 1, for typical values of the parameters taken here as 6± = 1, 
ai = 1.1 and /3±i = 10, the distance to threshold T (dashed line) given by Eq. (1461) and the 
non-dimensional nonlinear coupling coefficient A = A/(n 0 Tj_i) (solid line), with A given by 
Eq. (I571h as a function of 9\\. This graph is typical of the general behavior of these functions 
and shows that they are both decreasing as 9\\ increases, with A possibly reaching negative 
values, but only below threshold. In order to show that the value A c , given by Eq. (1581) . of 
A at threshold is positive in a wider range of parameters, we display in Fig. 2, as a function 
of /3±i for 9± = 0.2 (solid line), 9± = 1 (dotted line) and 9± = 5 (dashed line), the quantity 
min (A c ) obtained after minimizing A c in an interval of values of a p between 0 and a p i(/3_u). 
The latter quantity is arbitrarily defined such that the threshold is obtained for a value of 
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FIG. 1: Variation with of the distance to threshold T given by Eq. (|46l) (dashed line) and of 
the normalized nonlinear coupling coefficient A (solid line) evaluated from Eq. (1571) for 9± = 1 , 
a* = 1.1 and /3_|_j = 10. 



FIG. 2: Variation with /3_Li of the minimum min (A c ) of the normalized nonlinear coupling coefficient 
taken in an interval of values of a p between 0 and a p i(/3±i), defined such that the threshold is 
obtained for a value of equal to 100, for 6± = 0.2 (solid line), 6± = 1 (dotted line) and 8± = 5 
(dashed line). 

#H equal to 100. This graph shows that min(A c ) varies little with #j_ but is very sensitive 
to /3±i. As the latter parameter is increased, min (A c ) decreases towards zero but remains 
always positive. Although this numerical observation is definitively not a rigorous proof, 
it convincingly shows that A should remain positive in the parameter regime of physical 
interest. 

Thus, we can see that in the case of the bi-Maxwelian distribution functions for both 
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ions and electrons (i) the asymptotic model has the same structure as in the cold case, 
and (ii) it predicts the formation of magnetic holes which is defined by the sign of the 
coupling coefficient A . If the disributions are different for the bi-Maxwellian ones we can 
expect change of the sign of Aand appearance of magnetic structures in the form of humps 
respectively. In the next sections, we show how such mirror structures can be found for 
arbitrary distributions for both electrons and ions based on the variational principle when 
both pressures are functions of the magnetic held amplitude only. 

V. VARIATIONAL PRINCIPLE FOR STATIONARY ANISOTROPIC MHD 


As we saw in the previous sections, the nonlinearity for mirror modes originates from 
equations (l30j) which represent the anisotropic MHD in a static regime supplemented by 
corrections due to the FLR effects. Secondly, another origin of the nonlinearity comes from 
the drift kinetic equations, in particular, for the asymptotic model (1171) it comes from the 
stationary kinetic equations (1471) . Thus, the static anisotropic MHD together with the sta¬ 
tionary drift kinetic equations describe the nonlinear development of the mirror modes and 
its possible saturation in the form of static structures. In this section, we give formulation of 
the variational principle for such structures and establish connection it with the free energy 
formalism developed for the asymptotic model (1171) . 

A. Gyrotropic pressure balance 

We start from the pressure balance equation for a static gyrotropic MHD equilibrium 

0 = -V • P + i [j x B], (61) 

where the current j is defined from the Maxwell equation as j= -^VxB, and the pressure 
tensor P is assumed to be gyrotropic. The solvability conditions read B • (V • P)=0, and 
j-(V-P)=0. 

In terms of the tension tensor (6^ — bfij ) + S\\bibj, Eq. (IfiTh takes the divergence 

form £rn tJ = 0 where S± = p± + B 2 /(8 it) and Ay = p\\ — B 2 /(8tt), and the perpendicular 
and parallel pressures p± = T, a p± a and p± = T, a p± a are the sum of the contributions of 
the various particle species a. They are expressed as p± a = m a B 2 f pf a dv\\dp and p\\ a = 
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m a B f v^f Q dvjjd/u, in terms of the distribution functions f a , which satisfy the stationary drift 
kinetic equations 


V\\^\\fa 


h V ||^ + — V ||0 

'' “rv 


dfa 

dvw 


= 0. 


(62) 


(63) 


These equations are supplemented by the quasi-neutrality condition 

/ f a dv\\dn = 0 , 

a J 

that allows one to eliminate the electric potential. 

We consider partial solutions of the stationary kinetic equations (l62li which are expressed 
in terms of two integrals of motion: the energy of the particles W a = ujj/2 + fiB + (e a /m a )(f) 
and their magnetic moment p. In general, the solution can also depend on integrals which 


label the magnetic held lines 


34]. The choice f a = f a (W a ,p), as it will be shown further, 


can be matched with the solution found perturbatively for weakly nonlinear mirror modes 
within the asymptotic model (H7ll . In this case, the parallel and perpendicular pressures 
for the individual species and also the total pressures are functions of B only. We write 
P±a = P±a(B) and p\\ a = p\\ a (B). As seen in the next subsection, this property plays a very 
central role in the forthcoming analysis. 


B. Identity in the parallel direction and variational principle 


f- H, 


According to Ref. [33], the anisotropic pressure balance equation (jHTTl can be easily refor¬ 
mulated as follows: 


VP|| 


B 


(P± -P||)V£ = 


B x 


vx[?y=ia + 1 


Hence projection along the magnetic held gives 

4tt (p L-p\\) 


V lfP|| 


B 2 


B 2 


= °- 


47T 


B 


(64) 


(65) 


which coincides with Eq. (9.2) of Shafranov’s review [29]. It is possible to prove that Eq. 
(j65lh being solvability condition to (IbTTh reduces to an identity, by means of the statioinary 
kinetic equations (1621) together with the quasi-neutrality condition (1631) . To our knowledge, 
hrst time this_fact_was established by J.B. Taylor J 3 I, 32| and later by many others (see, 


for instance, 


33 


36]). Since the pressures depend on B only, Eq. 


reduces to 


dB 


(P±~P\\) 


B 


( 66 ) 
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In this partial case the system (1641) is written as 


B x 


V x 


(P± 


+ sH B 


= 0. (67) 

means that for station- 


V B 2 

The existence of the identity (RJ51) and its partial formulation 
ary states, the pressure balance provides only two scalar equations which, together with the 
condition V ■ B = 0, leads to a closed system of three equations for the three magnetic held 
components. 

It can be easilty shown that Eq. (1671) can be written in the following variational form: 


„ sb 

Bx — 

5 A 


= 0 


( 68 ) 


where T is given by the expression 

and A is the vector potential: B = [V x A] . In the pure 2D geometry with B = (B x , B y , 0) 
when the vector potential A has only one non-zeroth component ^ (^-component) for de¬ 
scription of stationary state we arrive at the variational principle hJ r =0, formulated in our 
paper 45). It is evident also that we have the same variational principle for stationary 
structures in r — z geometry when B = (B r , B z , 0). 


Note, that Eq. 


can be written also as 

P± ~ P\\ 


V x 1 + 47 t- 


B 2 


B 


= \B 


(69) 


where for scalar function y we have the equation (B ■ Vy) = 0. This equation shows that y 
is constant along each magnetic line. If the line is not closed so that at r —> oo B tends to 
the constant magnetic held B 0 and besides there both pressures p± and pyare also constant 
then for all such lines y = 0. Indeed, as we will show below, the equation for the stationary 
structures following from the guiding center formalism coincides with Eq. (1691) at y = 0. 


C. Derivation of the varitional principle from the guiding-center formalism 


To derive the variational principle for stationary mirror structures a three-dimensional, 
previously established for 2D configurations |45|, we now employ the Hamiltonian theory of 


guiding-center motion as stated in Section III of Ref. 


591. 
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Let us first consider a sort of particles with mass m and electric charge e. Instead of 
the particle position x and velocity v, we introduce new coordinates in the phase space: 
position X of the guiding center, parallel velocity component u along the magnetic held 
B(X, t), magnetic moment /i m\vj_\ 2 /2B(X, t), and a gyroangle (. The dynamics of 
these new unknown functions is determined by the following approximate Lagrangian (see 


derivation in Ref. 


59j). valid to the lowest order on spatial derivatives, 




- A(X, t) + mwb(X, t ) 
c 


• me • TYl 

x+ ~ t] - e0(Xj ^ (70) 


where A(X, t) is the vector electromagnetic potential, 0(X, t) is the scalar electric potential, 
b = B /B is the unit tangent vector. We see that Q is a cyclical variable in this (adiabatic) 
approximation, and therefore fi is a nearly conserved quantity. 

It will be important that the volume element in the non-canonical phase space (X, u, /i, £) 
contains a non-constant Jacobian J, 


me 


dV = dxdv = J(X^u)d i XdudndC t oc [B -|-w(b • cmlh^d^Xdud/idC,/(2n). 


Another important formula determines velocity of guiding center in stationary fields: 


• 7TICU C 

X ~ uh -|-—[curl b — b(b ■ curl b)]-wgrad UiB + ecp) x b. 

eB eB 


(71) 


It follows from Lagrangian (17U1) . This formula shows that the particle moves along the 
magnetic line (the first term), the second term is the drift velocity due to the centrifugal 
force (curlb — b(b ■ curlb) = [b x (b • V)b] where (b • V)b is the curvature); the last term 
is the drift due to the mirror force and the electric force. Eq. (1711) can be found in many 


papers, see, for example, 


34]- 


We consider (quasi-)stationary distributions of the given sort of particles like that 


dN = f (e(x, v), /i(x, v)) dV — [B -f- u(b • curl b)]F'(e, /i)d 3 Xdud/i<iC/(27r), (72) 

e 

where £ = /jB + e<p + mu 2 /2 is the Hamiltonian of a guiding center, and F(e,/j,) is a 
prescribed function of the two variables. It is assumed that F < 0 while F' e oc / > 0, and 
F —> 0 as £ —» +oo. It is clear that / satisfies the (collisionless) drift kinetic equation, since 
it depends on the exact integral of motion £ and on the approximate integral of motion /i 
(adiabatic invariant). Therefore there is no need in checking the hydrodynamic stationarity. 
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We require only two relations to close the model in a self-consistent manner: they are the 
Maxwell equation for a stationary magnetic field and the quasi-neutrality condition: 


X B jtotal/Cj 


(73) 


Ptotal 0, 

where j to tai and p to tai are the densities of the electric current and of the electric charge, 
respectively, produced by all sorts of particles present in the system. 

In the lowest order on gradients, the current density from the given sort of particles is (it 
follows from Eq.(3.53) of Ref.j59|) 


j/c = —V x (biV (p)) + (e/ c)iV(X). 


(74) 


Here N(p) = |M| = / pf Jdudpd(, where M is the spatial density of the magnetic moment. 
Using distribution (1721) . we have 

- b N(p) = -B J /iU(c, p)dpdu = -B^ J F(e, p)dpdu = B-^ (^j , (75) 

where p\\(B, <p) is the parallel pressure of the given sort of particles, 

p\\ (B,(f>) = B J mu 2 F' e (£, p)dpdu = —B J F(£,p)dpdu. 

It is remarkable that the calculation of N(X) = f X/ Jdudpd( with the help of Eqs. flTTT) 
and (1721) results in the following compact expression, 

(e/c)N(X) «Vx (bp,, /B). (76) 


Let us now label each 
substitution of Eqs. (175|) 
over a looks as follows, 


where p\\(B,<j>) = E« Pa 


sort of particle by an index a. Then the Maxwell equation (l73l) after 
and (1761) into Eq. (174l) for each a and after subsequent summation 


Vx { b B"^ P||(B ' 0) 

is the total parallel pressure, 


0, 


(77) 


p\\( B A) 


F (a){ £ a, p)dpdu, 

rv J 


with 


e a = pB + e a (p + m a u 2 / 2. 
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The quasi-neutrality condition 


J2 e aB / —F (a )(e a , n)dndu = 0 
is easily seen to have the form 

^P||(£,0) = O. (78) 

Since divB = 0, equations (|77l) and (17H1) possess the variational structure, and the corre¬ 
sponding functional is 

I = J[B 2 / (87t) — p\\(B, 0)]d 3 X. (79) 

In principle, the quasi-neutrality condition (17H1) allows one to express the electric potential 
(j) through B, and then the parallel pressure in Eq. f|79l) can be understood as a function of 
B only. As the result, we have the equation 


V x b 


B 


which coincides with Eq. 


variational principle previously derived in 


.47T 

= C 

). 

in 

45] 


~P\\{B) 


0, 


(80) 


at x — 0- Thus, we get a 3D generalization of the 2D 


45] by a different approach. 


It is worth noting that the quantity 

-^j =4vr[V x (bp\ { (B))} 

can be connected with mean (per volume unit) magnetic moment of plasma M=bpj|(i?), so 
that the magnetic field H = B + 47tM (in accordance with the definition of the Maxwell 
equations in continuous media. In this case, equation (I80p is nothing more as the Maxwell 
equation V x H = 0. 

Let us consider the most physically interesting case where the functions F( a ) (e a , //) have 
the exponential on e a form, 

F{a){e a ,ii) = - exp (-£ a /T a )D a (p), 

with constant temperature parameters T a and some positive functions D a {fi). In this case 
the u-integration is simple, and 

r- l-oo 

P\\{B,4>) = 5^)T a exp(-e Q 0/T Q ) / exp (-pB/T a )D a (p)dp, 
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where D a (/i) oc D a (p) by an a-dependent factor. Suppose we deal with the simplest electron- 
proton plasma. Then 

P\\(B,4>) = TiGi(B) exp(—e</>/Tj) + T e G e (B) exp(e0/T e ), 


where 


r -too 

G a (B) = B / exp (—/j,B/T a )D a (/j,)d/j,, a — i,e. 

Jo 

The quasi-neutrality condition (T78h now takes a simple form, 

exp (-e<j>/Tj)Gi(B) - exp (e<j>/T e )G e (B) = 0, 

from which we have (compare with Eq. ()53l) ) 

lii[G.(B)/G e (B)] 
p (i/Ti + i/jy ’ 

P\\(B) = (T t + T e )[G t (B)}-A[a e (B)]-A. 


(81) 


In particular, we may assume purely thermal isotropic electron velocity distribution, 
which corresponds to D e (/j,) = const. In that case G e (B) = const , and the total parallel 
pressure simplifies to 

G t {B) 


P\\ (B) = n 0 (Ti + T e 


Gi(B { 


o) 


Ti+T e 


(82) 


VI. TWO-DIMENSIONAL STATIONARY STRUCTURES OF THE GRAD- 
SHAFRANOV TYPE 

In two dimensions, we define the stream function 0 (or vector potential), such that 
B x = d'lp/dy, By = —d^/dx. In terms of 0 and B z , 

1 dB 2 z <90 


V x B 


x B 


= e 


( 


2 dx dx 


A 


!**)-’■{*'**>• (83) 
where {0 ,B Z } denotes the Jacobian. Furthermore, Vj_ = V--^Bi(Bi ■ V) — -§fe z (Bj_ • V), 
where V = (d x ,d y ) and B^ = (. B x ,B y ). 
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In Eq. (1671) , we now separate the (x, y)-components: 

-Vp_L + ' V )PJ- 


2 B 2 
1 

Air 


B 2 

(P±-P\\) 
An . 

1+ 32 ^ 


V - ^ B i(Bi • V) 


B 2 


(84) 


P\\) 


-VB 2 
2 2 


WipApEj = 0 . 


Due to identity 


the equation for the z component can be written 


B 
An 

1 

+ — 

47T 


4:71 

B ± V ) [l + ^(p±^p\\) 


i 4tt 
1 + ^(^ 


■Pll) 


(Bj_ • V)B~ — 0. 


In terms of ip, after integration, it leads to 


An 


1 + 


47T 

1b 


(p± - Pll)) = /WO- 


(85) 


( 86 ) 


Interestingly, in the isotropic case (p \ — V\\ = 0), we have B~ = B z (ip), in full agreement with 
the Grad-Shafranov reduction 126 


271 . 1291] . Furthermore, because the projection of the full 
equation on B is equal to zero, in the 2D case where the fields are functions of x and y only, 
the projection of Eq. (1841) on B^ vanishes identically. Therefore the relevant information is 
obtained by taking the vector product of Eq. (1841) with B^, in the form 


(87) 



B 2 1 

<1 

[ p± + 


) ~ (l Hr l ‘ v B ~ B) 


(ff 2 - B\ 
An 


An 

l + -^(p±-P||) 


A ip. 


This equation is supplemented by relation 

Equation (1H71) can be viewed as analogous to the Grad-Shafranov equation, the main 
difference being that the pressures are here prescribed as functions of the magnetic held 
amplitude. In particular, it does not reduce in the isotropic case to the usual Grad-Shafranov 
equation. Note that, according to the previous section, the obtained equations (1561) and (1571) 
follow from the variational principle for T. In particular, for the purely two-dimensional 


geometry when B z = 0 and B 2 = |Bj_| 2 Eq. (1571) reduces to 

V.{l + %( P± - Pll ) V*} = 0 (88) 

and thus derives from the variational principle 5J- = 0 with T = 4- f g(\Wip\ 2 )dxdy. Here 
the function g is found by integrating 

g\B 2 ) = 1 + ^{P± -P||)- 
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Due to identity 


, we have 
7 = 


87r 


P|l ) dxiy = - j Il^dxiy. 


(89) 


It follows that all the two-dimensional stationary states in anisotropic MHD are stationary 
points of the functional 7. Its density is a function of B = |V^| only. In the special case 
of cold electrons, this free energy turns out to identify with the Hamiltonian of the static 
problem 

Equations similar to (1HH1) arise in the context of pattern structures in thermal convection. 
As shown in [s 0 |, such equations represent integrable hydrodynamic systems. As in the usual 
one-dimensional gas dynamics, these systems display breaking phenomena where the solution 
looses its smoothness at finite distance, due to the formation of folds. As a consequence, 


these models require some regularization. For patterns, the authors of 


30] supplement in the 


equation an additional linear term involving a square Laplacian. In our case, this procedure 
corresponds to the replacement of 7 by 7 + (z//2) / (A ip) 2 dxdy, with a constant v > 0. In 
plasma physics, regularization can originate from finite Larmor radius (FLR) corrections, 
which are not retained in the present analysis based on the drift kinetic equation (see, e.g. 


XL, 


M5|): 


18|). In the three-dimensional geometry the same regularization reads as (compare to 

B 2 


7 = 


-P \\( B ) 


IV x B| 


7r. 


(90) 


_8vr J||V 7 2' 

One more remark. Let B be a function of x and y, but B z 7 ^ 0. In this case B z is not defined 
by stream function ip and therefore one needs to write down 



„ 6x11 

B x 

[ m 


= 0. 


(91) 


Hence it is easily to get Eqs. (9,10) from |45|. It is necessary to mention also that for the 
2D case the stationary states with B z 7 ^ 0 are determined from the equation 

67] 


V x 


hB 


For instance, the equation for B z has the form 

P±~P\\ 


1 + 47T- 


B 2 


= 0. 


B z = const, 


where instead of arbitrary function of ip (see Eq. (9) from 45]) we have const. In r — z 
geometry the analogous situation takes place where B v plays the same role as B z in the 
planar case. 
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Note that for isotropic plasma ( p± = pn), H. Grad and H. Rubin 28] formulated for the 
stationary MHD states the variational principle for 


7 = 


If - p ) *■ 


A. KP soliton 


We shall now show that the functional 7 we previously introduced has the meaning of a 
free energy. In the weakly nonlinear regime near the MI threshold, the temporal behavior 
of the mirror modes can be described by a 3D model 17, 18, Em]], that in the present 2D 
geometry reads 


u t = -\K 


I — 

5u 


with the free energy 


F = 


i(-£W 2 + u-^-u + (V±u) 2 ) + ^u 3 


dr. 


(92) 


(93) 


Here u denotes the dimensionless magnetic held fluctuations and e the distance from MI 
threshold. The third term in F originates from the FLR corrections, and A is a nonlinear 
coupling coefficient which is positive for bi-Maxwellian distributions. In Eq. (l92]h the 
operator | k y \ is a positive definite operator (in the Fourier representation it reduces to \k y \), 
so that Eq. (j92l) has a generalized gradient form. 

Let us now show that this result can be obtained from the functional 7 defined in (1891) . 
We isolate the perturbation <p in the stream function 0 = — Bq{x + p ) with ip —y 0 as 
|r| —> oo, so that the mean magnetic held B 0 is directed along the p-axis. We then expand 
Eq. (j89l) in series with respect to u. For the sake of simplicity, we restrict the analysis to 
the case of cold electrons. The expansion of the integrand B 2 /(8tt) — p\\ in F has then the 
form 


noT\\ 


(u T 1)^ 1 + u 

/3|l 1 + au_ 

= n 0 T|| [(/^i]" 1 - l) + u (a + 20I] -1 - l) 

+u 2 (—a 2 + a + — u 3 a 2 (a — 1) + ....] 


(94) 


where we use the usual notation /3m = 8tt7i 0 T\\/Bq. 
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As well known (see 


, e.g. 3, 


18]), near threshold, MI develops in quasi-transverse direc¬ 


tions relative to Bo- This means that, in the 2D geometry, p x ^ ip y and, with a good accu¬ 
racy, u coincides with < p x . However, in the expansion of u = \J(p x + l) 2 + p y — 1 — (p x +p y /2, 
it is necessary to keep the second term, quadratic with respect to ip. The linear term in the 
expansion of T vanishes and the quadratic terms is given by 


To = n 0 Tj| 


a(a-l) + ^ ]yl 


+ 


CL — 1 “f- 


2 1 




]} dxdy. 


where the factor a(a—l)+l//3 = —e/2 dehnes the MI threshold a = 1+1//3 _l (that the present 
equations of state accurately reproduces). It is also seen that for |e| -C 1, <p x /<p y ~ |£|~ 1//2 , 
in agreement with the quasi-one-dimensional development of MI near threshold. In this 
case, T 2 coincides with the quadratic term in (l93lh up to a simple rescaling and to the FLR 
contribution, Furthermore, the cubic term in (1941) gives the nonlinear coupling coefficient 
A = a (a — 1)>0. As a consequence, T , introduced in the previous section, reduces to the 
free energy of the asymptotic model. The temporal equation for ip has also the generalized 
gradient form originating from (192jh 

<Pt = ~ r T-with T =(95) 

op kz 


for which the associated stationary equation reads 


&P XX T Pxxxx ^Pyy Xd x (96) 

where the linear operator L = —ed xx + d yy — d xxxx is elliptic or hyperbolic depending on the 
sign of £. For £ > 0 (above threshold), this operator is hyperbolic, while below threshold it 
is elliptic and thus invertible in the class of functions vanishing at infinity. Remarkably, in 
the latter case, Eq. (19(1 identifies with the soliton for KP equation called lump. In standard 
notations, lump is indeed a solution of the stationary KP-II equation, 


Tu xx T u xxxx u yy T 3 {u ^ xx 0, 


(97) 


we see that 


where V is the lump velocity. When comparing this equation with 
plays the role of the lump velocity V and Xp x —y —3 u. 

The lump solution was first discovered numerically by Petviashvili 
now known as the Petviashvili scheme (see the next section). The analytical solution was 


57] using the method 
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later on obtained in 


58j . In onr notation, it reads 


12|e| (3 + e 2 y 2 — \e\x 2 ) 

A [3 + e 2 y 2 + |£|x 2 ] 2 

This function vanishes algebraically at infinity like r~ 2 . In the center region 
— \e\- 2 \e\x 2 — 3 < y < \e\~ 2 ^\£\x 2 — 3, the magnetic field displays a hole with a minimum 
at x — y — 0 equal to —4|e|/A. In the outer region, the magnetic lump has two symmetric 
humps with maximum values |e|/(2A) at y — 0 and x = ±3|e| —1//2 . The main contribution 
to the “skewness” I = / <~p x dx dy comes from the hole region, providing a negative value to 
/, in complete agreement with |17I . Il 8 ]. 


VII. NUMERICAL 2D SOLUTIONS 


In the 2D case, our regularized model equation for stationary pressure-balanced structures 
has a variational form 


(1 + ip x ) dg 

— d v 

1 

1_ 

(1 + u ) du 

y 

(1 + u ) du 


+ vA 2 ip 


0 . 


(98) 


Clearly, Eq. ([98]) describes stationary points 5J 7 /Sip = 0 of the functional T = J[g(u) + 
(z//2)(A ip) 2 ] dx dy , with some constant parameter v (in this expression and everywhere be¬ 
low, we use dimensionless variables). 

We applied two numerical methods to solve Eq. (1981) . The first one is a generalization 
of the well known gradient method which corresponds to a dissipative dynamics along an 
auxiliary time-like variable r of the form ip T = —T(5J r /S(p), with a positive definite linear 
operator T. It is clear that attractors in the phase space of the above dynamical system are 
stable solutions of Eq. (198 p . Unstable solutions however cannot be found by this method. 

Furthermore, the linear part of Eq. (1951) is of the form Lip = —g"(0)ip xx — g'(0)(p y y + vA 2 ip. 
The coefficient g"( 0) is proportional to £ (introduced in the previous section) and g'( 0) is 
positive within the adiabatic approximation. When these two are positive, the operator L 
is elliptic and it is possible to employ the so-called Petviashvili method [57]. It is a specific 
method for finding localized solutions of equations of the form Mip = N[ip], with a positively 
definite linear operator M and a nonlinear part N[ip], Note that in our case the Fourier 
image of M is 

M(k x , k y ) = g"(0)k 2 x + g'(0)k 2 + v[k 2 x + k 2 y ) 2 > 0. (99) 
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FIG. 3: Fig. 1. Unstable localized solution for v = 0.0004, r = 7, B s = 0.5 (in units of Bq), and 
C = 0.002. The value 1//3| = 1.127 prescribes an aspect ratio \Jg" (0)/ g' (0) = 0.2. 


In its simplest form, the iteration scheme of the Petviashvili method reads 

-7 


Vn+1 = (M IV [</?„,]) 


/ <p n Mip n dx dy 


( 100 ) 


f ip n N[(p n ] dx dy J 

where 7 is a positive parameter in the range 1 < 7 < 2. The corresponding multiplier 
strongly affects the structure of attractive regions in the phase space. 

It is worth noting that if the operator L is hyperbolic, solutions of the problem are not 
localized with res pec t to both x and y coordinates, and will be periodic or more generally 


quasiperiodic 


48 


50] • 


A. The results 

We performed computations with both numerical methods using fast Fourier transform 
numerical routines for the evaluation of the linear operators. Periodic boundary conditions 
for a computational square 2n x 27 t were assumed. 

For the gradient method, we used the simplest first-order Euler scheme for stepping along 
r, with 5r ~ 0.01. The operator T was taken in a form giving stable computation, namely 
F (k x , k y ) = l/[k 2 x + k 2 y + v(k 2 x + A;J) 2 ]. 

As for the Petviashvili method, the value 7 = 1.8 was used, leading, after an erratic 
transient, to a convergence of the iterations to unstable solutions of the variational equation 
§ 8 }. 
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FIG. 4: Fig. 2. Formation of a stable ID solution in a gradient computation, for the same 
parameters as in Fig.l. 

The main results of our computations can be formulated as follows. There do exist 
unstable localized solutions of Eq. (|98j) . which are similar to the lump solutions of KPII 
equation, when written in terms of u — d x ip (Fig. [3]). For asymptotically small e, they 
accurately coincide with KP solutions, independently of the electron temperature, as it 
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should be. Such low-amplitude stationary states do not depend on the particular choice 
of the regularization of g(u). No other kinds of solutions were found with the Petviashvili 
method. 

When the gradient method is used, large amplitudes u ~ 1 are achieved in many cases, 
and the final result turns out to be dependent on the choice of the parameters r and C in the 
regularized function g. Without regularization, no smooth stationary state is approached. 
Instead, a singularity occurs. Differently, when a regularized g with parameters r ~ 10 
and C ~ 0.001 is used, the final state identifies with a one-dimensional stripe in the form 
of a magnetic hole, as shown in Fig. []] that also displays typical stages of the “gradient” 
evolution. In all simulations, the magnetic held in the stripe was smaller than the ‘singular’ 
magnetic held B s given by Eq. (I56[) . For increasing r, the magnetic held in the stripe 
tends to decrease, down to 0. For initial conditions in the form of a slightly perturbed 2D 
lump, the hnal result is always a one-dimensional stripe of hole type, which demonstrates 


the instability of the 2D lump, in full agreement with the analytical prediction EH8|. 


In no cases stable 2D structures localized both in x and y directions were found. Instead, 
the gradient method showed that stable structures can only be one-dimensional, transverse 
to the magnetic held. An initial localized perturbation of sufficiently high amplitude develops 
into an increasingly long structure along the y axis, and eventually reaches the boundary of 
the computational domain. 

The question arises whether the ID shock solutions obtained in jgj (for which mini? > B s ) 
would identify with the present solution when v —> 0, a limit which is unreachable in the 
present numerics. It is possible that the presence of the bi-Laplacian regularization leads 
to overshooting in the shock solution, resulting in the convergence towards solutions where 
mini? < B s . 


B. 2D mirror structures with B z / 0: stripes and magnetic bubbles 

Let us consider some numerical examples. For simplicity we take the function D^g) oc 
exp (gB s /Ti) for g < g*, and Dj{g) = const for g > g *, with some parameter B s < B 0 , and 
a large /i*. Such constant-like behaviour of Di(g) at very large g is necessary both from 
formal and physical points of view (see discussion in Ref. 45|). At B = B 0 we thus have 
a nearly Gaussian ion perpendicular velocity distribution with the temperature Tj_(i?o) = 
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FIG. 5: Some plots corresponding to expression (11011) . 


Ti /(1 — B s /Bo ) > T % . The distribution becomes strongly non-Gaussian as the magnetic held 
decreases to values B < B s . Let us normalize all magnetic held values to B 0 so that formally 
B 0 = 1. As the result, we have the following expression for the ratio Gi(B)/Gi(l) = W(B), 

B{ 1 - C)(l - B s ) {1 - exp [—/?(£> - B a )]} 


W{B ) = 


{1 - exp[—i?(l - B a ))} 


(B - B s 


+ C'exp[l?(l — B)], (101) 


with a sufficiently large regularizing parameter R and a small parameter C. Some plots, 
with B s = 0.4, R = 7.0, for several C, are shown in Figj5] 

We substituted this dependence into Eq. (l82l) and then into Eq. (j90]) . with T e = 0. To hnd 
stable stationary 2D mirror structures with B z ^ 0, we parametrized magnetic held in the 
following manner, 


B x = -il> y (x, y ), B y = y ), B s = 7(x, y). 


We hxed mean values (B x ) = 0, (B y ) = cos0n, (B~) = sin0 o . Then we employed the 
gradient numerical method described in Ref. 45J [with a simple generalization to include 
y(x, y)\ to hnd minimum of the functional 


Bon — 


RW + 7 2 ) + ^(|Ai/>| 2 + |V T | 2 )]d 2 X, 


( 102 ) 
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FIG. 6: Some examples corresponding to expression (11031) . Here B s = 0.4, R = 7.0, C = 0.003. 
System with 1//3 m = 0.2 is linearly unstable. System with 1//3|| = 1.7 is linearly stable, but 
subcritical structures are possible. System with 1//3|| = 5.0 is stable, no structures are possible. 


where 


g(B) = W(B), (103) 

P II 

and /3|| = 87rnoTi/B^. Plots of function g(B), for several values of /3y, are shown in FigjGj The 


mirror instability takes place when the second derivative g" { 1) is negative. Subcritical mirror 


structures are possible when g" (1) is positive, but there is a range of B where g "(B ) < 0. 

It is important that besides purely ID stable configuration (“stripes”), in our computations 
we have detected for some parameters also essentially 2D stable solutions — “bubbles”, as 
shown in Fig{7]for B s = 0.4, R = 7.0, C = 0.003, (B y )/B 0 = 0.2, l/f3\\ = 1.71. In general, 
“bubbles” takes place when B z dominates, i.e. cos@o is sufficiently small. They have the 
perfect circular shape in the case when B x = 0 and B y = 0 (see Figj8]). In all cases we 
have inequalities g” (B in ) > 0, and g" (B out ) > 0, so the unstable range of B is passed in the 
vicinity of the bubble boundary. When B± = 0 the magnetic fields are constant inside and 
outside circle everywhere accept transient layer which is defined by the FLR. The size of the 
circular patch is defined by two factors: the conservation of magnetic field flux and the cell 
size. The FLR introduces small input in the this constraint, it plays a role of the surface 
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tension. 

In Fig[9]is shown for circular bubbles the diagram of all possible both stable and unstable 
states at the fixed (3 n measured by the B s field. Because of the magnetic fields outside and 
inside the bubbles are constant, stability and instability of each state is defined by the second 
derivative of the function g(B). At the given (3\\ the Bi and B 2 curves represent the inner 
and outer magnetic fields when FLR is absent. The FLR in this case provides a transient 
solution matching the inner and outer regions. But to say that these are the inner or outer 
solution one needs to have another jump, or some patch if we speak about two-dimensional 
structures. Both states Bi and B 2 are linearly stable. These states satisfy the necessary 
boundary conditions, namely, continuity of the magnetic field: H\ — H 2 = H , where H is 
an additional constant. These states, thus, can be considered as conjugated states, or, by 
another words, these are bistable states. When changing (3\\ which is defined has a meaning 
of the parameter e we move along the curves Bi and B 2 . One should mention that in this 
case /3|| is some auxiliary dimensionless parameter. Real (3\\ is found depending on a state 
by means of B\ or B 2 . If one considers any state B, say, at the given f3 \\, without any 
conjugation, then one can get linear stability or linear instability by analyzing the second 
derivative sign of the function g(B). The second point is that by fixing two conjugated 
states one can say only that B 2 > B\. Only in the case when you have another jump one 
can say whether it is a hole or a hump. One more point is that the case considered here 
corresponds to the pure B z case when B± = 0. 

VIII. CONCLUSION 

In the first part of this paper we presented a review of our results concerning the weakly 
nonlinear regime of the mirror instability in the framework of the so-called asymptotic 
model. This model was demonstrated to belong to the class of the gradient systems for 
which the free energy can decrease in time only. In particular, it was shown that the sta¬ 
tionary localized solutions of the model, below the mirror instability, occur unstable and, 
above the threshold, the system has a blow-up behavior up to amplitude comparable with a 
mean magnetic field that is typical for subcritical bifurcation. We showed also that account 
of electrons (increase their temperature) does not change the structure of the asymptotic 
model. For bi-Maxwellian distribution functions for both electrons and ions all analyzed 
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FIG. 7: Example of 2D “bubble”, with 1//3|| = 1.71, {B v )/Bq = cos0q = 0.2. 
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FIG. 8: Circular “bubble”, with 1//3|| = 1.71, cos0q = 0.0. 



FIG. 9: B 1 -B 2 diagram of stationary states depending on /3m measured by B s 
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structures within the model have the form of magnetic holes. Humps can appear for dis¬ 
tributions different from the bi-Maxwellian ones. For instance, such situation is possible 
after a stage of quasi-linear relaxation ( for details see results of numerics 3.7]). The second 
part of this paper contains original results concerning the possible two-dimensional mirror 
structures which can be formed at the saturation regime of subcritical bifurcation. In par¬ 
ticular, a detailed analysis was presented for the Grad-Shafranov equations describing static 
force-balanced mirror structures with anisotropic pressures given by equations of state de¬ 
rived from drift kinetic equations, when assuming an adiabatic evolution from bi-Maxwellian 
initial conditions. It turns out that in two dimensions, the problem is amenable to a varia¬ 
tional formulation with a free energy provided by the space integral of the parallel tension. 
Slightly below the mirror instability threshold, small amplitude solutions associated to KPII 
lumps are obtained and shown to be unstable. Based on the variational computation (the 
gradient method) of the stationary mirror structures, this instability is shown to result in 
appearance of one-dimensional stripes when the magnetic fields outside and inside stripes 
are homogeneous with a jump which structure is defined by the FLR effects. Such two- 
dimensional evolution of the stationary structures are formed for below and above threshold 
of the mirror instability when the R 2 -component of magnetic field is absent. For the finite 
but small enough values of B z the resulting structures represent stripes. With increasing B z 
instead of stripes we observed in numerical simulations the formation of magnetic bubbles 
with the homogeneous magnetic field inside the bubbles. When B z becomes larger B± the 
form of bubbles change their form from elliptic to the circular one when B± = 0. In the 
latter case, the magnetic field outside and inside bubbles occurs constant and undergoes 
jump due to the FLR effects while crossing the bubble. In this case, the FLR effects play 
the role of surface tension. Note also, when considering stable subcritical structures, the 
drift kinetic approximation breaks down, as the deep magnetic holes obtained by a gradient 
method appear to be strongly sensitive to the regularization process, an effect which in a 
more realistic description could be provided by FLR corrections and/or particle trapping. 
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